function centers = gaussian_fit_3d(image, sigma1, sigma2, threshold)
    % 读取3D图像并转换为double类型
    image = double(image);
    
    % 计算高斯差分（DoG）
    gauss1 = imgaussfilt3(image, sigma1);
    gauss2 = imgaussfilt3(image, sigma2);
    dog = gauss1 - gauss2;
    
    % 归一化
    dog = dog / max(abs(dog(:)));
    
    % 检测局部极大值
    local_max = imregionalmax(dog);
    peaks = (dog > threshold) & local_max;
    
    % 获取极大值点的位置
    [rows, cols, slices] = ind2sub(size(peaks), find(peaks));
    centers = zeros(length(rows), 3);
    
    % 进行3D高斯拟合
    for i = 1:length(rows)
        x = cols(i);
        y = rows(i);
        z = slices(i);
        
        % 提取局部区域进行拟合
        window_size = 5;
        [X, Y, Z] = meshgrid(x-window_size:x+window_size, ...
                             y-window_size:y+window_size, ...
                             z-window_size:z+window_size);
        
        sub_volume = dog(max(1, y-window_size):min(size(dog,1), y+window_size), ...
                         max(1, x-window_size):min(size(dog,2), x+window_size), ...
                         max(1, z-window_size):min(size(dog,3), z+window_size));
        
        % 构造拟合数据
        X = X(:); Y = Y(:); Z = Z(:);
        sub_volume = sub_volume(:);
        
        % 3D 高斯拟合模型
        gauss_model = @(b, X) b(1) * exp(-((X(:,1)-b(2)).^2/(2*b(3)^2) + ...
                                          (X(:,2)-b(4)).^2/(2*b(5)^2) + ...
                                          (X(:,3)-b(6)).^2/(2*b(7)^2)));
        
        % 初始参数 [幅值, x0, sigma_x, y0, sigma_y, z0, sigma_z]
        init_params = [max(sub_volume), x, 1, y, 1, z, 1];
        
        % 拟合优化
        options = optimset('Display', 'off');
        fitted_params = lsqcurvefit(gauss_model, init_params, [X, Y, Z], sub_volume, [], [], options);
        
        % 记录中心点位置
        centers(i, :) = fitted_params([2, 4, 6]);
    end
    
end
